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We investigate minimal energy solutions with vortices for an interacting Bose-Einstein condensate 
in a rotating trap. The atoms are strongly confined along the axis of rotation z, leading to an effective 
2D situation in the x — y plane. We first use a simple numerical algorithm converging to local minima 
of energy. Inspired by the numerical results we present a variational Ansatz in the regime where 
the interaction energy per particle is stronger than the quantum of vibration in the harmonic trap 
in the x — y plane, the so-called Thomas-Fermi regime. This Ansatz allows an easy calculation 
of the energy of the vortices as function of the rotation frequency of the trap; it gives a physical 
understanding of the stabilisation of vortices by rotation of the trap and of the spatial arrangement 
of vortex cores. We also present analytical results concerning the possibility of detecting vortices by 
a time-of-flight measurement or by interference effects. In the final section we give numerical results 
f^ . for a 3D configuration. 
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After the achievement of Bose-Einstein condensates in trapped atomic gases M many properties 
of these systems have been studied experimentally and theoretically |2| . However a striking feature 

■^l- . of superfluid helium, quantized vortices p), p), has not yet been observed in trapped atomic gases. 

■^j- ' There is an abundant literature on vortices in helium II, an overview is given in M. 

The atomic gases have interesting properties which justify efforts to generate vortices in these 

\& • systems: the core size of the vortices is adjustable, as in contrast to helium the strength of the 

interaction can be adjusted through the density; the number of vortices in atomic gases can be in 
principle well controlled; for a small number of particles in the gas metastability of the vortices 
can be studied, that is one can watch spontaneous transitions between configurations with different 
number of vortices. 

S Several ways to create vortices in atomic gases have been suggested. A method inspired from liquid 

helium consists in rotating the trap confining the atoms S] ; at a large enough rotation frequency it 
becomes energetically favorable at low temperatures to produce vortices; two different paths could 
be in principle followed: (1) producing first a condensate then rotating the trap, or (2) cooling the 
gas directly in a rotating trap. It has been recently proposed in p] to use quantum topological effects 
to obtain a vortex. Other methods that do not rely on thermal equilibrium have been suggested M , 



Here we study theoretically the minimal ener gy configurations of vortices in a rotating trap M . 
The model is defined in section O; in sections III to VI we assume a strong confinement of the 
Vh ' atoms along the rotation axis z so that we face an effective 2D problem in the transverse plane 

x — y. We present numerical results for solutions with vortices that are local minima of the Gross- 
Pitaevskii energy functional (section |HJ). These solutions contain only vortices with a charge ±1, 
the vortices with a charge larger than or equal to 2 are thermodynamically unstable (section |fV|). 
We discuss possibilities to get experimental evidence of vortices in atomic gases in section M. Finally, 
we concentrate on the regime where the interaction energy is much larger than the trap frequencies 
uj x ,y, the so-called Thomas-Fermi limit M. This is complementary to the work of [hfj. We obtain in 
this "strong interacting" regime analytical pred ictions base d on a variational Ansatz that reproduce 
satisfactorily the numerical results (section (VJ) . In section VII we present results for vortices in 3D, 



that is in a trap with a weak confinement along the rotation axis. 



II. MODEL CONSIDERED IN THIS PAPER 

The atoms are trapped in a potential rotating at angular velocity Q. In the laboratory frame the 
Hamiltonian of the gas is therefore time dependent. To eliminate this time dependence we introduce 
a rotating frame at the angular velocity Q so that the trapping potential becomes time independent; 
this change of frame is achieved by the single-atom unitary transform: 



U{t)=e n ' Lt/h (1) 

where L is the angular momentum operator of a single atom. As the unitary transform is time 
dependent the Hamiltonian in the rotating frame contains an extra inertial term, given for each 
atom by 

ihU 1 (t)j-U(t) = -Q ■ L. (2) 

The atoms are interacting via the effective low energy potential commonly considered in the 
literature, V{f\ — fa) = <73_D<5(ri —fa), where the coupling constant is related to the s-wave scattering 
length a (taken here to be positive) and to the atomic mass m by 330 = A-Kh 2 a/m. In this paper we 
consider the case of a dilute gas (with a density much smaller than a~' i ) at zero temperature. We 
can then assume that the N particles of the gas are condensed in the same state 4>. The wavefunction 
4>(f) is time independent as we are in the rotating frame and minimizes the energy per particle in 
the condensate, given by the Ginzburg-Landau energy functional: 

w]=y<^ m +- 2 ^r- ^ 

In this energy functional TCq contains the kinetic energy and the trapping potential energy of the 
particles: 

h 2 

«o = -;r-A + l7(?,t = 0). (4) 

2m 

The energy functional includes also the inertial term — Cl-L and a term proportional to |0| 4 describing 
the interaction energy between the particles in the mean field approximation. 

As done in the present experiments with atomic gases we take the case of a harmonic trap, with 
eigenfrequencies u a : 

U(f,t = 0)= ^2 2 mClJ ° r «- ( 5 ) 

a=x,y,z 



Furthermore in all but in section VII 



we will assume that the trapping potential is much stronger 
along the z axis than along the x, y axis, with an oscillation frequency much larger than the typical 
interaction energy Ng3D\4>\ 2 P er particle. This situation, although not realized experimentally yet, 
is not out of reach, in particular when one uses optical traps rather than magnetic traps |0|. In 
this strong confining regime the motion of the particles along z is frozen in the ground state of the 
strong harmonic potential: 

<j>(x,y,z) ~ ii{x,y) (— r^-J e"™* 2 /2S (6) 

and only the dependence of the wavefunction ip in the x — y plane remains to be determined. 
Inserting Eq.ffl) in the energy functional Eq.(H) we get the corresponding energy functional for ip to 
be minimized, dropping the constant term ^hui z : 

E ^^ = J dr m + w (r) 

This 2D energy functional gives the energy per particle in the condensate measured from the zero- 
point energy along z. The 2D Hamiltonian is 



j. = -— A x>y + - ^2 mwlrl. (8) 
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The trap is now rotated around the z axis at the angular velocity Q. The interaction term 
involves an effective 2D coupling constant between the atoms [H2j 



9 = 93D {^h) ■ (9) 

Most of the results of the paper are dealing with the 2D energy functional; a nu merical result for a 



We concentrate on 



local minimum of the full 3D energy functional will be given in the section VII . 

the so-called Thomas-Fermi regime, where the interaction energy per particle is much larger than 

hbJx,y The opposite regime has already been studied in J10|. 



III. LOCAL MINIMA OF ENERGY WITH VORTICES 

In this section we briefly discuss the general problem of minimizing energy functionals of the type 
Eq.(M). We present the numerical algorithm that we have used and we give numerical results for 
the 2D problem. 



A. A numerical algorithm to find local minima 

The algorithm in our numerical calculations is commonly used in the literature to minize energy 
functionals E[ip,ip*] of the form Eq.(|7|). The intuitive idea is to start from a random ip and move 
it opposite to the local gradient of E[ip, ip*] that is along the local downhill slope of the energy. 
Numerically this is implemented by an evolution of ip parametrized by a fictitious time r: 

-^ =*>•*']• < 10 > 

Assuming a ip normalized to unity we get the following equation of motion for ip: 

-^-iP = [h±-ql z + n 9 \iP\ 2 - I i{t)]iP (li) 

CLT 

that is a non-linear Schrodinger equation in complex time t — —ir. The quantity fj, appearing in 
this equation can be expressed explicitly in terms of a functional of ip, ip*; it ensures that the norm 
of ip does not evolve with r. 

This equation is, for fi = 0, standardly solved by a splitting technique, propagating during the 
time step dr first with potential energy in position space (where it is diagonal) then with kinetic 
energy in momentum space (where the Laplacian is diagonal). One goes back and forth between 
position and momentum space with Fast Fourier Transforms along x and y. In our case fi ^= and 
the Hamiltonian contains L z = xp y — yp x ', we have therefore complemented the splitting scheme by 
(1) a propagation during dr due to — TiQ.xp y in position space along x and momentum space along 
p y , and (2) a similar procedure for the hQ,yp x propagation, that is in momentum space along x and 
position space along y. 

One can check that the mean energy of ip is a decreasing function of r: 



, n , , , - 

—E[iP,r] = ~2 d 2 r 



5E 



SiP* 



<0. (12) 



In the case we consider E remains always positive: the scattering length and therefore g are positive 
so that the interaction term is positive, and |fi| is smaller than the trap frequencies uj XtV so that the 
centrifugal potential — imfi 2 r 2 cannot exceed the trapping potential. Therefore E has to converge 
to a finite value for r — » oo. Asymptotically dE/dr = and ip satisfies j^[ip,ip*] — 0, so that we 
recover for r = oo the time independent Gross-Pitaevskii equation: 

fiiP = T-C±ip + Ng\ip\ 2 iP - nL z iP (13) 

where fi — /x(r = oo) is now the chemical potential of the gas B. 

As we have started from a random wavefunction ip, without assuming any symmetry properties of 
ip, we expect the trajectory ip(r) to converge as r — + oo to a local minimum of the energy functional. 
We have checked this assumption by adding a small random wavefunction to ip and resuming the 
evolution in r; ip was relaxing to its initial value. Mathematically the steady state solutions for the 
t evolution that we find are stable, which is equivalent to saying that they are local minima of the 
energy. Note that not all solutions of the Gross-Pitaevskii equation share this property: the Gross- 
Pitaevskii equation expresses only the fact that the energy functional is stationary in ip, which is 
the case e. g. at saddle points of the energy functional (an example is a vortex with a charge \q\ > 1, 



see section IV) 



B. Numerical results in 2D 

Applying the algorithm detailed in the previous section we present results on local minima of the 
energy functional Eq. Q7J) for asymmetric and symmetric traps in the x — y plane. We characterize 
the non-axisymmetry of the trap by e such that 



uj x =uj/{l + e) (14) 

w„=w(l + e). (15) 

In FigJll we show different local minima configurations obtained for e = 0.3 and a rotation frequency 
fl = 0.2u>; each configuration has been obtained for different random initial ip's. The holes observed 
in the spatial density correspond to the vortex cores. We have always found that the phase of ip 
changes by 2tt around a vortex core; we have not found vortices with a charge ±q, where the integer 
q is strictly larger than one; this fact will be explained in the next section. Furthermore the sense 
of circulation is the same for all vortices. 

To quantify the effect of the non-axisymmetry of the trap we have plotted in Fig.H the dependence 
of energy of different vortex configurations on u) x /u} y for a fixed ui; we measure the energies from 
Ei ao , the energy of the zero- vortex solution in the axisymmetric case e = 0. The zero- vortex solution 
exhibits a significant variation of energy with e; for a non-zero e the wavefunction ip develops a phase 



proportional to Q, for weak D.'s, which accounts for the energy change as explained in section VI B , 
The solutions with vortices experience quasi the same energy shift as function of e. As only the 
energy difference between the various local minima matters we will from now on only consider the 
axisymmetric case e = to identify the solution with the absolute minimal energy. 

Note that the solutions ip with several vortices obtained in the limiting case e = are not eigen- 
vectors of L z ; this reflects a general property of non-linear equations such as the Gross-Pitaevskii 
equations to have symmetry broken solutions; it is explained in [hoi how to reconcile this symmetry 
breaking with the fact that eigenvectors of the full iV-atom Hamiltonian are of well defined angular 
momentum. 



IV. STABILITY PROPERTIES OF VORTICES 

In this section we recall that a (normalized) wavefunction tp such that E[ip,ip*] has a local min- 
imum in ip, describes a condensate having all the desired properties of stability, that is dynamical 
and thermodynamical stability. We then show that a vortex centered at r — with an angular mo- 
mentum strictly larger than ft is not a local minimum of energy and is therefore thermodynamically 
unstable. 



A. Stability properties of local minima 

Let us express the fact that ip corresponds to a local minimum of the energy. A first condition is 
that the energy functional is stationary for ip, that is ip solves the Gross-Pitaevskii equation Eq.(|l3|). 
To get the second condition, we consider a small variation of ip, 

iP^iP + Sip (16) 

preserving the normalization of the condensate wavefunction to unity: 

\\ip + 5ip\\ 2 - HVII 2 = = (ip\Sip) + (5ip\ip) + (Sip\6ip). (17) 

We expand the energy functional E[tp, ip*] in powers of Sip, neglecting terms of order Sip 3 or higher. 
Using Eq.(|l7|) and Eq.(|l3|) we find that terms linear in Sip vanish so that 

6E = i«*# WIKc ( jJjW \ +o{8rp 2 ). (18) 

We have introduced the operator 

Hgp + Ng\iP\ 2 Ngip 2 

NgiP* 2 H* GP +Ng\i>f 

and the Gross-Pitaevskii Hamiltonian 

Hgp = H± + Ng\i[,\ 2 - Q,L Z - At. (20) 

The fact that E has a local minimum in ip imposes that the Hermitian operator C c be positive. In 
general C c will be strictly positive apart from the zero energy mode {iip,—iip*) corresponding to 
an inessential change of the global phase of ip. We now show that the positivity of C c implies the 
stability of the solution ip. 



C <= ( »r„.J ^ Tr.,.,,2 ) ( 19 ) 



1. Dynamical stability 

Consider first the problem of so-called "dynamical stability" : to be a physically acceptable con- 
densate wavefunction, ip has to be a stable solution of the time dependent Gross-Pitaevskii equation 

iW t i\> = Hgp^ (21) 

otherwise any small perturbation of ip, e.g. the effect of quantum fluctuations or experimental noise, 
may lead to an evolution of tp far from its initial value. To determine the evolution of a small 
deviation Stp as in Eq.(|lfj|) we linearize Eq,(]2l|): 

«*(KMi!$) (22) 

where the operator £ is related to C c by 

As ip is time independent, so is C and dynamical stability is equivalent to the requirement that the 
eigenvalues of £ have all a negative or vanishing imaginary part. As we now show the positivity of 
C c leads to a purely real spectrum for C. Consider an eigenvector (u, v) of C with the eigenvalue e. 
Contracting Eq.(|23|) between the ket (\u), \v}) and the bra ((u|, (v\) we get 

e[(u\u)-{v\v)} = ((u\,{v\)£ c (MY (24) 

Note that the matrix element of C c is real positive as C c is a positive hermitian operator. We now 
face two possible cases for the real quantity (u\u) — (v\v): 

(u\u) — (v\v) — 0. In this case £ c has a vanishing expectation value in (\u), \v)); as C c is positive (|w), \v)) has to be an 

eigenvector of C c with the eigenvalue zero; from Eq.([23|) and the fact that I is invertible we find that (\u), \v)) 

is also an eigenvalue of C with the eigenvalue 0, so that e — is a real number. 

(u\u) — (v\v) > 0: we get e as the ratio of two real numbers, so that e is real. 

2. Thermodynamical stability 

A second criterion of stability is the so-called "thermodynamical" stability. For zero temperature, 
this condition can be formulated in the Bogoliubov approach B, where the particles out of the 
condensate, which always exist because of the interactions, are described by a set of uncoupled 
harmonic oscillators with frequencies esign [(u|w) — (v\v)] /%, where (u, v) is an eigenvector of C with 
the eigenvalue e. In order for a thermal equilibrium to exist for these oscillators, their frequencies 
should be strictly positive, which is the case here in virtue of Eq.(E4J) [h3[. If a mode with a negative 
frequency were present thermalization by collisions would transfer particles from the condensate tp 
to this mode, leading to a possible evolution of the system far from the initial state if) JlJ] . 

What happens for solutions ip of the Gross-Pitaevskii equations that are not local minima of 
energy ? The operator L c has at least an eigenvector with a strictly negative eigenvalue. In this 
case one cannot have thermodynamical stability, that is one cannot have e \{u\u) — (v\v)] > for 
all modes Jl3(. From the non-positivity of C c one cannot however distinguish between a simple 
thermodynamically instability or a more dramatic dynamical instability. 

B. Why not a vortex of angular momentum larger than h ? 

For simplicity we consider only a single vortex in the center of an axi-symmetric trap. We 
show that vortices with a change of phase of 2q-K are not local minima of energy, that is are (at least 
thermodynamically) unstable. We have found numerically a solution of the Gross-Pitaevskii equation 
Eq.(n3) by an evolution in complex time, starting from a wavefunction ip with an angular momentum 
qh along z, as already done in 113]; our solution of the Gross-Pitaevskii equation with imposed 
symmetry is a local mimimum of energy in the subspace of functions with angular momentum qh 



along z, but not necessarily a local minimum in the whole functional space, as we will see for 
\q\ > 1. In the Thomas-Fermi regime p S> hui we find that the solutions can be well reproduced by 
a variational Ansatz of the form 



i>(x,y) 



iqO 



tanh K q r]' q ' 



/' - 



1 2 2\ 1/2 



Ng 



where 8 is the polar angle in the x ~ y plane and where jl, the chemical potential in the lab frame 

jl — fi + qhQ. 

does not depend on Q. In this Ansatz the vortex core is accounted for by tanh' 9 ', a function that 
vanishes as r' 9 ' in zero as it should, and the condensate density outside the core coincides with 
the Thomas- Fermi approximation commonly used for the zero- vortex solution M. We calculate the 
mean energy Eq. Q7I) of the variational Ansatz and we minimize it with respect to the variational 
parameter K q \ we get 



fim 



1/12 



where 



du u (tanh 



*M 



(u)-l)' 



is a number (ci = 0.7687, c 2 = 0.5349, . . .). 

In order for the vortex of charge q to be a local minimum of energy, the operator L c of Eq.(|19|) 
has to be positive. This implies that the operator on the first line, first column of £ c , the so-called 
Hartree-Fock Hamiltonian, be positive: 

Hhf = H±+ 2Ng\iP\ 2 -fr + qMl- QL Z > 0. 

To show that this is not the case it is sufficient to find a wavefunction f(x, y) leading to a negative 
expectation value for TLhf- As the potential appearing in Hhf has a dip at r — we have taken / 
of a form localized around r — 0: 



(25) 



(26) 



(27) 



(28) 



(29) 



f(p,v) = 



cosh 



n^ 



1/2 



where 7 is adjusted to minimize the expectation value. For e.g. q — 2 we take 7=1 leading to 

(f\H HF \f) 



(f\f) 



-0.407/2 + 2M1 



As Q < uj <C [i. this quantity is negative. A similar conclusion is obtained for q > 2. 

We have also performed a numerical experiment, evolving tp in complex time starting from Eq. 
we find that the vortex q = 2 splits in two vortices q = +1 symmetrically dispatched |l4| . A nu- 
merical diagonalization of C shows that the vortex q = 2 is alternatively dynamically and thermo- 
dynamically unstable when one increases jl/hco Jig], 



(30) 



(31) 



V. HOW TO DETECT THE PRESENCE OF VORTICES ? 



Several signatures of the presence of vortices have been proposed in the literature. A first possi- 
bility is a measurement of the excitation spectrum as studied in [I17J. Another idea is to measure 
the second order correlation function of the atomic field [[18| . 

A third signature of the presence of vortices is also the holes in the density due to the vortex 
cores. As the size of the vortex core in the Thomas- Fermi regime is too small to be observed in 
situ by optical imaging techniques, we suggest to switch off the trapping potential and let the cloud 
expand; as we now check the size of the cloud and the size of the vortex cores are magnified by the 
same factor in the expansion, so that the cores become observable. 

To study the expansion of the gas when the trap in the x — y plane is switched off, the confinement 
along z being kept constant, one has to solve a 2D time dependent Gross-Pitaevskii equation. In this 



section the trap is axi-symmetry with a time dependent frequency ui(t). We consider the evolution 
in the laboratory frame, as the detection is performed in this frame: 



ihdtipiab — 



-— A + -mcj 2 (t)r 2 + Ng\i{>i ab \ 2 
2m 2 



tplab- (32) 



As shown in ]19| , P0| the effect of the time dependence of u>(t) can be absorbed by a scaling and gauge 
transform of the wavefunction: 

V^(r, t) = j^e mr2 A/2fi ^(r/A(t), t) (33) 

where ui(0) is the oscillation frequency before opening the trap; the scaling parameter solves: 

A=^- w a (t)A (34) 

with initial conditions A(0) = 1, A(0) = 0; if the trap in the x — y plane is abruptly switched off at 
t = + the scaling parameter is given by 

A(i) = y/l+w 2 (0)t 2 . (35) 

Introducing the renormalized time r given by 

■Wi = dT (36) 

we find that tp solves the same equation as ipi a b with a constant trap frequency equal to u>(0): 

h 2 . l 



ihd T tp ■■ 



-^-A + >w 2 (0)r 2 + Ng\i>\ 2 
2m 2 



i). (37) 



As ipiab rotates in the trap at the frequency SI in the lab frame, so does ip in terms of the renormalized 
time t. In the limit of t — ► oo, r tends to a finite value r max , so that ip is rotated by a finite angle 
during the ballistic expansion: 

QT -* = n J A% = I^y- (38) 

Therefore %pi a b rotates with respect to its value when the trap is switched off and its size is magnified 
by \{t). 

A fourth possibility, giving direct access to the phase of the vortex, is to measure the interference 
fringes between a condensate with vorticity and a reference condensate with no vortex. We study 
this possibility as an application of the scaling solution [El| . The condensate 2 has one or several 
vortices and is originally centered at f = 0, the condensate 1 has no vortex and is centered initially 
at r = d. After ballistic expansion of the condensates the resulting density can be written: 

Pl+2 = \ p \/\irn(r-3)V2 M + p l/2 e ,^/2« e , S ( r1e , 7 (t) |2 (3g) 

where pi and p-2, are the densities of the condensates 1 and 2 respectively, S is the phase due to 
vorticity of the condensate 2, y(i) is a relative phase depending only on time; to obtain the phase 
terms quadratic in r we have used Eq.(B3|) with the asymptotic value A/A ~ 1/t for t — > oo. We 
have plotted an example of interference fringes with two vortices in Fig.H. 

The above scaling result is exact only for axi-symmetry 2D traps. For a non axi-symmetric 
traps and for 3D situations where the confinement along z is not strong, it has been shown that 
approximate scaling solutions exist in the absence of a vortex (ly,M ; in presence of vortices we have 
integrated numerically the time dependent Gross-Pitaevskii equation and found that the density 
experiences an approximate scaling, the magnification of the vortex being slightly larger than the 
one of the cloud (see also Q). 



VI. INTUITIVE VARIATIONAL CALCULATION 



To get a better understanding of the numerical results we now proceed to an intuitive Ansatz 
for the wavefunction with several vortices. It coincides very well with the numerical results and 
allows an easy construction of the minimal energy configurations with vortices. It gives a physical 
understanding of the stability conditions and of the structure of the solutions: a set of n vortices is 
equivalent to a gas of interacting particles in presence of an external potential adjusted by the rota- 
tion frequency of the trap. We restrict to the cas e of an axi-symmetric trap, a good approximation 
for weak (< 10%) non-axisymmetries (see section IIIB). 



A. Ansatz for the density 



To construct the Ansatz we split ip in a modulus and a phase: 



ip(x,y) = \tp\e . 

In the Thomas-Fermi regime, the modulus in presence of n vortices appears as a slowly varying 
envelope given by the Thomas-Fermi approximation used in the 0-vortex case: 



<Pslow 



1 2 2 

/i — ^muj r 



Ng 



1/2 



with narrow holes digged by the vortices with charge q — ±1, represented by tanh functions of 
adjustable widths and with zeros at adjustable positions: 



(40) 



(41) 



\tp\ = ipsiow x TTtanh[«fc|f- a k R\]. 



(42) 



The positions of the vortex cores dk are expressed in units of the Thomas-Fermi radius R of the 
condensate: 



R 



2// 



From section 



IV B 



we expect as typical values for the inverse width of the vortex cores Kk ~ 
(nifi/h 2 ) 1 ' 2 . The chemical potential is not an independent variable but is expressed as a function 
of the other parameters from the normalization condition {^\il>) — 1; neglecting overlap integrals 
between the holes we get 



M = Mo 



l + 2]T(l-aI) 7 i^ + 0( 7 -L l! 



(ACfei?) 2 



-(kR) 4 



where 1/(kR) 4 ~ (Tuo/fi) 4 <C 1 and where fj,o is the Thomas-Fermi approximation for the condensate 
chemical potential without vortices: 



Ho = 



mu> 2 Ng 



1/2 



(43) 



(44) 



(45) 



B. The phase 

The general form of the phase of tp in Eq.(Ud) in presence of n vortices is: 

n 

S(x,y) =^qk8k +So(x,y) 



(46) 



where the integer q^ = ±1 is the vortex charge (that is the angular momentum (over h) of the vortex 
k with respect to its core axis), 9t is the polar angle of a system of Cartesian coordinates (X,Y) 



centered on the vortex core and So is the single-valued part of the phase. The function So can in 
principle be determined from the modulus of ip from the continuity equation: 

Aw{\ifj\ 2 v] = 0. (47) 

The local velocity field v is related to the phase S by 

v=— VS-fiAr. (48) 

m 

This expression is derived from the relation between the velocity operator and the momentum 
operator in the rotating frame, v = p/m — O A f . Expanding the continuity equation we obtain 

\4>\ 2 AS + V|V| 2 ■ VS - ^[xd y - yd x M 2 = 0. (49) 

This can be turned into an equation for the single- valued part So of the phase; because the density 
\ip\ 2 in a trap vanishes at the border of the condensate So is uniquely determined (up to a constant) 
by the resulting equation (see Appendix); this is to be contrasted to the case of superfluid helium 
in a container, where the flux, not the density, vanishes at the border, which requires a boundary 
condition on the gradient of the phase. 

Eq.([f3) can be solved for a non-axisymmetric trap in the absence of vortices. The solution is 
given by 

n u x + U)y 

which leads to a change in the energy per particle 

5E = --rrrtf ,) l '*>- (51) 



6 ' ' (w| + Wy)WJW 



y 



where [itf is the Thomas-Fermi approximation for the chemical potential for Q — 0, ^tf — 
(mcOxLOyNg/iv) 1 ' 2 |23). As can be seen in Fig.PJ this prediction is in good agreement with our 
numerical results. 

In presence of vortices the equation for S is more difficult to solve analytically. From now on we 
consider the case of an axi-symmetric trap, as the energ y ord ering of the vortices solutions is not 



affected for weak (< 10%) non-axisymmetries (see section [II B). For a single vortex at the center of 



the trap one can see that So = solves Eq.(W!l). From the spatial dependence of the phase obtained 



numerically (section IIIB) for a displaced vortex or several vortices we have identified the following 
heuristic Ansatz, obtained in setting cu x = to y — uj in Eq.(BOj): 

S o (x,y) = (52) 

that we will use in the remaining part of the section. 

C. Further approximations for the mean energy 

In the calculation of the mean energy, we make some further approximations in the spirit of the 
Ansatz Eq.(|42J). The reader not interested by these more technical considerations can proceed to 
the next subsection. 

The kinetic energy involves an integral of the gradient squared of the wavefunction: 

|VV>| 2 = W 2 [(Vln |t/>|) 2 + (VS) 2 ] . (53) 

For the gradient of the modulus of tp we neglect the variation of the slow envelope ip 3 i ow : 

n . 

Vln|^|^^K,^[ Kfc |r-a fc i?|](e r ) fe (54) 

fc=i 

where (e r ), — (r— (XkR)/\r — dkR\- The terms in this sum are peaked around the vortex cores; 
assuming a separation between the vortex cores much larger than their width, we neglect all the 
crossed terms in the square of Eq.(p4|). Consider now the second term in Eq.(B3f). The gradient 
squared of S involves diagonal terms (W^) 2 and non-diagonal terms V6k-V6 k / ; the modulus squared 



of ip involves holes with a density varying as 1 — tanh 2 = sech 2 . In the following we keep the sech 2 
for the vortex k only if it is multiplied by (V9k) 2 , a quantity diverging in the center of the core; the 
other terms lead to converging integrals smaller by a factor (/i/hco) 2 , which is the inverse surface of 
a vortex core ( f d 3 r \ip a i ow \ 2 sech 2 Kr oc 1/(kR) 2 ). This finally leads to 



h 2 

Ekin ^ g — / d 2 r \i) s i ow \ 2 



2_]tanh 2 [K fc |r - a k R\](V8 k ) 2 + k\ (tanh'[K fc |r - a k R\\)' 



In the same spirit we simplify the contribution of —QL Z to the energy: 



E rot = -Kl ■ I d 2 r U\ 2 rA VS 



—HQ- I dr\if} s i ow \ FA } y v</ ; . 



E™* 



In the potential energy 

E pot — / d 2 f I' 
we will neglect in \ip\ products of sech coming from different vortex cores. 



-muj 2 r 2 + -Ng\ip\ 2 



(55) 

(56) 
(57) 

(58) 



D. A more physical form of the mean energy 

After the approximations detailed in the previous subsection the mean energy in presence of n 
vortices is a sum of one-vortex self energies and binary interaction energies between the vortices: 



n 

E = -mo + E H/ ( Q "fe' Kfe ) + o 



2 EE 

fc = l k'^k 



V(dk,dk' 



Eq.(p9|) allows to interpret a system with n vortices as a gas of "particles" with binary interactions; 
the form of the interaction potential obtained here is valid for a separation between the "particles" 
larger than the size 1/k of the vortex cores. The vortex self-energy is 



W(a, k) 



{!+<»" 



(M 2 

Mo 
-qhQ, [(1-a 2 ) 2 ] 



C - In kR ■ 



1 



ln(l 



Mo 
3 



(41n2- 1 



(l-« 2 ) 2 
{nRY 



where C = 0.495063. The lines in Eq.(pCJ) correspond successively to Ekin, E ro t and E pot . This can 
be seen as an effective potential for the vortices. One can check that the part of W independent 
of Q expells the vortex core from the trap center, whereas the part proportional to fi provides a 
confinement of the vortex core (see the following subsection). 
The vortex interaction potential is given by 



1 - Ti 

-V(a,J3) = —q a qf3 \ d 2 r \^ alow \ 2 V6 SR ■ V0g 



iiR ■ 



This interaction term is equivalent to the one found in the homogeneous case and describes a repulsive 
interaction for vortices turning in the same direction (q a q0 > ) and is attractive for vortices 
with opposite charges M. An attractive interaction will lead to the coalescence and consequently 
annihilation of vortices with opposite charges. Therefore we find in stationary systems always 
vortices with equal charges. 

As the interaction potential V(a, j3) does not depend on the parameters k we can optimize sepa- 
rately the self-energy part with respect to k and find 



(59) 



(60) 



(61) 
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where £ — [§(4 hi 2 — 1)1 ~ 1.08707. By rewriting the above equation as 



ftV 



\e 



(i — —muj {aR,y 



we find that k 2 is proportional to the local chemical potential at the position aR of the vortex core. 
We finally get the explicit form for the self energy as 



W(a) 



(hu) 
Mo 



{i + (l- 2 ) 



2 In 2+ 1 vno , 2 , Q^o 

^3— + ln ^7 + ln(1 - Q) -^ 



2 (l-« 2 ] 



where v = 0.49312. 



E. Case of a single vortex: critical frequencies 



(62) 



(63) 



(64) 



In Fig.0 we have plotted the self-energy of a vortex as a function of the displacement of the core 
from the trap center, for different values of the rotation frequency Q. The analytical prediction 
coincides very well with the numerical value p4j . 

For $1 = the position of the vortex at the trap center gives an energy maximum. For Q > the 
rotation of the trap provides an effective confinement of the vortex core at the center of the trap for 
positive charges q (see the term proportional to Q, in Eq.(p4J)); from now one we therefore take all 
the charges q^ to be equal to +1. For a large enough f2 we reach a situation where a vortex at the 
trap center corresponds to a local energy minimum, by further increasing fi the vortex state at the 
trap center becomes a global minimum with energy less than the condensate without vortex. 

The above suggests that we have to distinguish two critical rotation frequencies: The first one 
defines the frequency Q s tab above which the vortex is a local minimum of energy. Above the frequency 
Q c the single vortex solution has an energy lower than the condensate without vortex. We calculate 
Qstab from the condition d 2 W/da 2 — at a = and Q c from the condition W — at a = 0: 



Q c — In 

Mo 



ft.. 



2/i 



hi 



Tiuj 



where C' — e 



(21n2+l)/3+l/2 



1.8011. As we are in the regime 



Qstab J25|. Our prediction for Q, c scale as (log/xo)//io as in 
better agreement with the numerics. 



o S> Tiuj Q c is approximately twice 
with a coefficient C' leading to 



(65) 
(66) 



By integrating Eq.(f 
with equal charges: 



F. Case of several vortices 



11) we get an explicit form for the vortex interaction potential for vortices 



V{a, (3) 



(hu 



fio 



a 2 + f3 2 



+ 2^ 



■ a ■ 0) log 



\a A /3|atan 
1 



-A/3| 



(l-a-f3) 
2a-j3 + u 2 p 2 ~ 



|q-/3| 4 



At short distances between the two vortex cores the logarithmic term in the above expression 
dominates, leading to a repulsive potential ~ —2(1 — a ■ p) log [a — [3\(hu)) 2 /fj,o- In Fig.ra we plot the 
interaction energy between a vortex at the center of the trap and one of equal charge displaced by 
aR; the interaction is purely repulsive. A conclusion which essentially holds as well for arbitrary 
vortex positions. In FigJrl we show the total (interaction + self-energy) for two vortices symmetrically 
displaced from the trap center, as function of the displacement; the analytical prediction coincides 
again very well with the numerical results [M . To obtain the equilibrium distance between the two 
vortex cores one minimizes the total energy over a in Fig.H. 



(67) 
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To get the minimal energy configurations as function of the rotation frequency of the trap, we 
minimize our analytical prediction for the energy over the positions of the n — 1, 2, . . . vortex cores. 
The result is shown in Fig.W. Each curve corresponds to a fixed value of n; it starts at O = Q sta b(n) 
(for Q < Qstab(n) there is no local minima of energy with n vortices); it becomes the global energy 
minimum for fi = fi c (n). We have plotted these two critical frequencies as function of n in FigJq. 

We have also given numerical results (circles) in Fig.M. Even if there is good agreement between 
analytical and numerical results, we still need a numerical calculation to check the stability of the 
solutions; our simple analytical Ansatz is indeed not sufficient to predict the destabilization of a 
given vortex configuration at high Q, a phenomenon studied with a numerical calculation of the 
Bogoliubov spectrum for a single vortex in ]2q ]. 

For a fixed value of the number of vortices n there may exist local minima of energy, in addition 
to the global minimum plotted in FigM a situation known from superfluid helium 0. E.g. for n — 6 
(see Fig.H) the global minimum of energy is given by a configuration with six vortex cores on a 
circle; there exists also a local minimum of energy with one vortex core at the center of the trap and 
five vortex cores on a circle. The energy difference per particle between the two configurations is 
very small, 5E ~ 0.002?kj for the parameters of the figure and probably beyond the accuracy of our 
variational Ansatz. For relatively large rotation frequencies Q one can find local minima of energy 
configurations with many vortices (see W for superfluid helium) ; we plot two configurations with 18 
vortices in Fig.hoL with an energy difference 5E = 0.0034?kj. 

In estimating the physical relevance of these energy differences one should keep in mind that 
NSE matters, rather than SE, where N is the number of particles in the condensate: e.g. at a 
finite temperature T the ground energy configuration is statistically favored as compared to the 
metastable one when NSE ~^> ksT. 

VII. VORTICES IN A 3D CONFIGURATION 

We have extended the numerical calculation to the case of a 3D cigar-shaped trap, that is with 
a confinement weaker along the rotation axis than in the x — y plane. Even in this case rotation 
of the trap can stabilize the vortex. We show in Figjlll density cuts of a solution with 5 vortices; 
the vortex cores are almost straight lines in the considered Thomas- Fermi regime, except at vicinity 
of the borders of the condensate. As in section IVll the core diameter is determined by the local 
chemical potential in the gas. 

This suggests that our 2D Ansatz (section |Vl| ) can be generalized to 3D situations, with a*, and 
Kk depending on z. 

VIII. CONCLUSION AND PERSPECTIVES 

We have presented in this paper an efficient numerical algorithm and a heuristic variational 
Ansatz to determine the local minima energy configurations for a Bose-Einstein condensate strongly 
confined along z and subject to a rotating harmonic trap in the x — y plane. 

Our results can be used as a first step towards finite temperature calculations. Interesting prob- 
lems are e.g. the critical temperature for the vortex formation and the Magnus forces induced by 
the non-condensed particles on the vortex core [ B7J . 
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APPENDIX A: UNIQUENESS OF THE PHASE FROM THE CONTINUITY EQUATION IN A TRAP 

Consider two solutions Si and S3 of the continuity equation: 

div[H 2 VS] = If-Qcdy - yd x )W 2 - (Al) 
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Si and S2 correspond to the same positions of the vortex cores, so that their difference S12 is a 
single- valued function of the position, solving 

div[|V>| 2 ^Si2] = 0. (A2) 

To show that in the case of a trapped condensate S12 is a constant we consider the following integral 

1= f f div[|V| 2 S12VS12] (A3) 



where the integration runs over the area A of the condensate. 

First, we transform / using Gauss's formula into an integral over the border A of the condensate: 



/fSiaVSia ■ n = (A4) 

J A 

which vanishes as l^l 2 = on the border of the condensate. 
Second, we expand the integrand of I as 

divO^SiaVSia] = SiadivOi/fVSia] + |V| 2 (VSi 2 ) 2 . (AS) 

The first term in the right hand side vanishes in virtue of the continuity equation. Therefore 

= /= (J M 2 (VS 12 ) 2 . (A6) 

As the integrand is positive this implies VS12 = 0, that is Si 2 =constant. 
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FIG. 1. Isocontours of the density \ip\ 2 for a 1-vortex (a) and a 4-vortex (b) configuration obtained numerically in a 
non-axisymmetric trap (e = 0.3) with a rotation frequency fl — 0.2a;, and /i ~ 40hui (see text for the definition of e, ui); 
the unit of length for x and y is (h/mtu) 1 ' 2 . 

FIG. 2. Numerical results for the energy of 0- vortex, 1-vortex and 4-vortex configurations (from top to bottom), as function 
of the frequency ratio ui x /uj y , for a fixed product of the frequencies u> 2 = lj x lo v . The rotation frequency of the trap is Q, — 0.2a; 
and the energies are measured in units of Tioj from the 0-vortex energy Ei SO in the axi-symmetric trap; the chemical potential 
is /i ~ AOTiuj. The circle s correspond to numerical results; the solid line is an analytical prediction for the 0-vortex case as 
obtained in section 



VI B 



FIG. 3. Isocontours of the total density pi+2 for two ballistically expanded condensates. The condensate 2 has two vortices; 
it was prepared in an axi-symmetric trap with Q — 0.2o>. The condensate 1 has no vortex and is slightly displaced along the 
axis x connecting the two vortex cores at time t. 

FIG. 4. Self-energy of a vortex in an axisymmetric trap as function of the distance aR of the core from the trap center, for 
[io ~ 80hu). (a) SI = 0.03a; and (b) Q, — 0.045a;. The solid lines are given by the analytical prediction W(a). The stars are 
obtained numerically. The critical frequencies defined in the text are il c = 2£l sta b — 0.06a;. Ei 30 is the energy of the 0-vortex 
solution and the unit of energy is hui. 

FIG. 5. Interaction energy between a vortex at the center of the trap and a vortex at a distance aR from the center, as 
given by the analytical formula Eq.(pj7|) for V. The trap is axisymmetric; the unit of energy is (Tito) 2 / no- 

FIG. 6. Energy of a system of two vortices symmetrically displaced by ±aR from the trap center, as function of the 
displacement a, for fl = 0.1a; and no — 80hu>. Solid line: analytical result. Stars: numerics. The trap is axisymmetric; Ei ao is 
the energy of the 0-vortex solution and the unit of energy is Tiui. 

FIG. 7. Local minima of energy with n vortices in an axisymmetric trap as function of the rotation frequency fl of the 
trap in units of the trap frequency u). The chemical potential no is approximately 40hu>. Note that for a fixed n we kept only 
the local minimum with the lowest energy. The circles are numerical results. The reference of energy is Ei 30 , the energy per 
particle in the absence of vortex, and the unit of energy is Tiuj. The lines with increasing absolute value of the slope correspond 
to n = 1, . . . , 4 vortices respectively. 

FIG. 8. Critical frequencies flstab(n) and f2 c (n) (see text) obtained from the analytical Ansatz as function of the number 
of vortices. The chemical potential no is approximately 40hu and the unit for Q is the trap frequency u>. 

FIG. 9. In an axi-symmetric trap with Q = 0.3a; and no — 40hcu, different configurations of 6 vortices corresponding to a 
local minimum of energy: (a) 6 vortices on a circle, with an energy per particle E = Ei SO — 0.5910fto;. (b) one vortex at the 
center and 5 vortices on a circle, with E — Ei SO — 0.5890?ia;. Ei ao is the energy per particle in the absence of vortex. The unit 
of length is (h/muj) 1 ' 2 . 
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FIG. 10. In an axi-symmetric trap with Q = 0.5w and /j, ~ AOhui, different configurations of 18 vortices corresponding to 
a local mininiuni of energy: (a) with one vortex at the center E = Ei SO — 2.3988hu>, this is the global minimum of energy, 
(b) without vortex core at the center; the energy per particle is slightly higher, E = Ei SO — 2.3954?icj. Ei ao is the energy per 
particle in the absence of vortex. The unit of length is (h/mco) 1 ' 2 . 

FIG. 11. A local minimum energy solution in 3D with 5 vortices obtained from the numerical evolution in complex time. 
The trapping frequencies are in the ratio u> x — lo v — 4u) z . The chemical potential is 53.7hu) x ,y The trap is rotated at a 
frequency Q — 0.25ui x ,y (a) Isocontours for a cut of the density in the plane z — 0, showing the 5 vortex cores, (b) Isocontours 
for a cut of the density in the x — z plane, showing the dependence with z of the vortex cores. The unit of length is (S/nw,,,) 1 ' 2 . 
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